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Abstract 

In this work, we propose a numerical method based on high degree continuous nodal ele- 
ments for the Cahn-Hilliard evolution. The use of the p-version of the finite element method 
proves to be very efficient and favorably compares with other existing strategies (C^ elements, 
adaptive mesh refinement, multigrid resolution, etc). Beyond the classical benchmarks, a nu- 
merical study has been carried out to investigate the influence of a polynomial approximation 
of the logarithmic free energy and the bifurcations near the first eigenvalue of the Laplace 
operator. 

Introduction 

We consider an isothermal binary alloy of two species A and B, and denote by w e [—1, 1] the 
ratio between the two components. By thermodynamic arguments, and under a mass conservation 
property, Cahn and Hilliard described a fourth-order model for the evolution of an isotropic system 
of nonuniform composition or density. They introduced a free energy density / to define a chemical 
potential, and use it in the classical transport equation (see [10], [12] and [l-^]). The total free energy 
of the binary alloy is a volume integral on Q of this free energy density (bulk free energy) : 

J" := I /(u,Vu,V2m, ...) dK (0.1) 

They assumed / to be a function of u and its spatial derivatives. A truncated Taylor expansion of 
/ has thus the following general form: 

f{u) - f{u) + L-\Iu + Ki® V^it + Vu-K2- Vw, (0.2) 

where V is the Nabla operator. By symmetry arguments, they showed that L = and Ki and K2 
are homothetic operators. Moreover they used Neumann boundary condition to cancel the term 
in V^w which yields 

T := j {f{u) + K\Vu\^)dV, (0.3) 

where k is a parameter (often denoted e^/2) which is referred to as the gradient coefficient. 
Then, the chemical potential w is defined by: 

w := f'{u) - 2kAu. (0.4) 

A is the Laplace operator. If we denote by J the flux and by Ai{u) the mobility, the classical Fick 
law provide the following equations: 

dtu^ -W ■ J a,nd J ^ -M{u)\/w. (0.5) 
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Finally, the Cahn-Hilliard equation takes the following general form: 

' dtu^V ■[Miu)Vw], on SI CM'', 
w = ip{u) - e^Au, on il C W^, 

Vu • = = Vui • I', on dft, 



(0.6) 



where t denotes the time variable, e (= V2k) is a measure of the interfacial thickness, ip (= /') is a 
nonlinear term, A4 is the mobility function, v is the outward pointing unit normal on the boundary 
dil. It is well known that the Cahn-Hilliard equation is a gradient flow in H^^ with Lyapunov 
energy functional J-. 

For a regular uniform alloy, the free energy / is explicitly given by: 
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where fcs is the Boltzmann constant, Nm a molecular density, T the temperature and Tc > T the 
critical temperature. Thus the nonlinear term ip is: 



f :ui-^ -NrnksTcU + In 

2 \1 — u 



(0.8) 



which is singular at li = ±1. These singularities give rise to the flrst difficulty in a numerical study, 
so this function tp is often replaced by the derivative of the classic quartic double-well potential, 
where / takes the following form: 

f -.u^^il-u'Y' , (0.9) 



with derivative: 



ip : u u 



(0.10) 



The Cahn-Hilliard equation has been extensively studied in the case where tp is replaced by a 
polynomial fimction (see [12], [22] and [31J]). Furthermore, this model has been used successfully for 
describing phase separation phenomena, see for example the survey [28[, and the references therein, 
or other recent results on spinodal decomposition and nucleation in [5, 7, 17, 25, 26, 32, 33, 37]. 
Recently, Ma and Wang have studied the stationary solutions of the Cahn-Hiliard equation (see 
123]). The case of non smooth ip has been the object of much less research (see JS] and [15[). 

Other frequent simplifications are often made. The mobility Ai is often assumed to be con- 
stant and the physical parameters are set to 1 - as we have done above in (0.9). For a more 
physically relevant choice of mobility, we mention [:>(>[ where the following form is proposed 
A4{u) = max{0, 1 — u^}. Among the physical parameters, e has a peculiar role since it may lead to 
different asymptotic behaviors and equilibria (see [24] and section 3). The study of evolution with 
e — >■ is of great importance: in particular a constant mobility leads to a Mulhns-Sekerka evolution 
(nonlocal coupling) whereas a degenerate mobility leads to a purely local geometric motion (see 
[4]). Furthermore, when the interface thickness is of the order of a nanometer, an artificially large 
parameter e is often used to regularize the numerical problem. When a fine resolution is out of 
reach, a change in the height of the barrier between wells in the free energy density, coupled with 
a change on e, allows simulations with larger length scales (see [35] for details). 

The evolution of the solution of (0.6) can essentially be split into two stages. The first one is 
the spinodal decomposition described in section 2 where the two species quickly separate from each 
other. In longer time, the evolution is slower, and the solution tends to reduce its interfacial energy. 
These two evolutions require different methods for an efficient global simulation. In the beginning, 
a very small time step and a precise grid resolution allow efficient computation. But this is not 
appropriate to get long-time behaviors. So an adaptative time accurate or/and an adaptative mesh 
can improve the efficiency of the algorithms. However, in the long-time evolution, the interfaces 
have to be precisely captured so that a global adaptative mesh cannot be used. In the literature, 
many technical ideas have been studied: adaptive refinement of the time-stepping or of the mesh. 
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elements (see [.)■")]), multigrid resolution (see [21]). 

We propose here an alternative method using high degree lagrangian nodal finite elements 
under a constant mobility A4 = 1. The use of p- version (increasing polynomial degree, see [2]) 
instead of /i-version (decrease mesh-step) has proved to be efficient for propagation [1, 19, 20], 
corner singularities ],^)4], or oscillating problems ]!)]. The numerical results obtained here with the 
finite element library Melina ]27] show that this method is suitable in the Calm-Hilliard framework 
as well. 

Our paper is organized as follows: in section 1, we shortly describe the discretization (in both 
time and space) including the nonlinear solver and the high degree finite elements we used. Section 
2 and section 3 are respectively devoted to the numerical results for the one-dimensional and 
the two dimensional problem. We investigate the performance of our method through different 
quantitative and qualitative aspects of the Cahn-Hilliard equation: comparison to explicit profile- 
solution in ID (see section 2), spinodal decomposition (see section 2), discussion about polynomial 
approximations of the logarithmic potential (see section 2). impact of the temperature and the 
parameter e (see section 2 and 3), long-time behavior and asymptotic stable states (see section 3). 
The numerical results are compared with existing ones in the literature, validating our approach. 



1 Discretization 

1.1 Space-Time schemes 

We start with the description of the time discretisation. Given a large integer N, a time step t, and 
an initial data {wq,uq), we denote by {wn,Un)n<N the sequence of approximations at uniformly 
spaced times t„ = nr. The backward Euler scheme is given by: 

= AWn+1, 



(1.1) 

Wn+1 = 1p{Un+l) - £^Aw„+i. 

A Crank-Nicolson scheme could easily be implemented but our experiences show that it gives results 
quite similar to the ones we shall show in the sequel. The schemes are immediately generalized 
to our case. We denote by (•, •) the scalar product in L'^{Q). We use the standard Sobolev space 
H^(ri) equipped with the seminorm 

\h\i^\\Vhh2, 

and with the norm 

I1^!1i = (|/^I? + I1/^IIlO'^'- 

The weak form of the equation (1.1) reads: 

{un+i - Un,x) = -t(A^(u„+i)Vw„+i, Vx), for all x e Xi, 

(1.2) 

{wn+1,0 = (V'(wn+i),0 + (£^Vu„+i, VO, for aU C e ^2, 

where Xi and X2 are the spaces of test functions (H^(ri) for example). We discretise in space by 
continuous finite elements. Given a polygonal domain fl, for a small parameter /i > 0, we partition 
Q into a set T'' of disjoint open elements K such that h = max (diam(X)) and II A' = 51. 

Thus, we define the finite element space 

- {x e C{n) -.xlKeP for aU K e r'^} , (1.3) 

where P is a space of polynomial functions, see section 1.3. We denote by the standard 

basis of nodal functions. Thus, for u and v € C(ri), we define the lumped scalar product by: 



(u,u)'' ^{u,ip,){v,ipj){ip,,ipj). 
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The scheme (1.2) can be rewritten m the fuUy discrete form, just by replacing the continuous scalar 
product with the lumped scalar product. 

We denote u = {uj)ji^j and w = {wj)j^j, the finite dimensional representation of u and if 
(we omit here the subscript n of the time scheme). Then we define the matrices A and M. whose 
coefficients are given by the following relations: 

[A]ij {WipijS/ipj), "stiffness" matrix, foralli,j€ J, 

[M]^^- := {ipi,(pj), "mass" matrix, for all i,j E J. 

For each time-step, given a previous solution (w„, u„), (w„-|_i, u„-|_i) is solution of the system 
rAw„+i +Mu„+i = Mu„, 

(1.4) 

Mw„+i -e^Au^+i -M*(u„+i) = 0, 

where ^ is a pointwise operator {related to tjj), and with (wo, Uo) the finite dimensional representa- 
tion of the initial data. The system (1.4) is clearly block-symmetric. The proof of the convergence 
of this scheme can be found in [3]. 



1.2 Nonlinear solver 

At each time step, we use a Newton procedure to solve the implicit nonlinear system (1.4) . For 
(1.4), we define the operator L by: 

tA M 



M 





Then denote by S the matrix of the left hand side of the backward Euler scheme, 

S = 

Denote also by G the following operator: 

-M*(u) 

Finally denote by Y„ the couple (w„,u„) for each n < N. The backward Euler scheme at each 
time-step satisfies the following formula: 

LY„+i + G(Y„+i) - SY„ = 0. (1.5) 

The Newton iterates (Y,^ := (w,'j, ufj)fegN satisfy for each n < N 

Y° - Y 

(1.6) 

Y^i = Yf, -[L + Dg (YfO] [(L + G - S) (Y,';)] , for all k e N, 

where Dq, (Y^) is the differential of G at point Y^. Actually, we stop the procedure at A: = fc„ 
when the residual is small, and define Y„+i := Y^". System (1.6) is an implicit linear system for 
each Newton-step, handled with a biconjugate gradient method. 

When the nonlinear term is logarithmic, we should deal with the singularities at ±1. However, 
in all our computations, the solution stays far from ±1 so that no special care is needed. This is 
expected. Indeed, it is known that in the one dimensional case the solution satisfies an L°° bound 
which is strictly less than one (see [14]). The same result has not been proved in higher dimension 
but it is probably true. 

A simple remark shows the mass conservation through the total scheme. Indeed, if we multiply 
the first component of the second equation in the system (1.6) by the vector I := (1, 1, 1) which 
belongs to V'^, we get for all fc £ N: 

IMu^ = IMu„. 
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1.3 Implementation with high degree finite elements 

The finite element library Melina [27] has the feature of providing lagrangian nodal elements with 
order up to 64 (the nodes may be chosen as the Gauss-Lobatto points to avoid Runge phenomenon 
for large degrees) . It can thus be used as a p- version code - see [2] - or even to implement spectral 
methods - see [(>]. In the following results, we use quadrangular elements for two-dimensional 
computations, with degree from 1 to 10. So we use the notation Qi with i £ {1, 2, 3, 4, 5, 6, 7, 8, 9, 10} 
to describe these elements. We justify this strategy by the fact that the expected solution is smooth 
but may present a thin interface ; since high degree polynomials are able to capture high frequencies 
they are well suited in such situations. Some comparisons are shown below between degree 1 on 
a refined mesh, and degree 10 on a coarse mesh, justifying the efficiency of the method (in both 
terms of accuracy and computational cost). 

2 Cahn-Hilliard evolution. Polynomial approximation of the 
logarithm 

The temperature plays a crucial role in the evolution of the solution. The function ip defined in 
(0.8) depends on two values of the temperature T and Tc- When the temperature T is greater than 
the critical temperature Tc, the second derivative of ip is non-negative, thus function -0 is convex 
and has only one minimum. We say that the function ip has a single well profile. Thus the solution 
tends to this unique minimum and the alloy exists in a single homogeneous state. 

But when the temperature T of the alloy is lowered under the critical temperature Tc, the 
function ip changes from a single well into a double well (see Figure 1), and the solution rapidly 
separates into two phases of nearly homogeneous concentration. This phenomenon is referred to as 
spinodal decomposition. If the initial concentration belongs to the region where the energy density 
is concave, i.e. between the two spinodal points (t_ and a-i- (see Figure 1), the homogeneous state 
becomes unstable. 

The concentrations of the two regions composing the mixture after a short stabilization have 
value near the so called hinodal points /3_ and j3j^ (see also Figure 1), defined by 

/'(/?_) = /'(/3+) = /(/^+)-/(/^-) ^ ^ith /?_ < /?+. (2.1) 

P+ — P- 

If the free energy is symmetric, the binodal points are the minima of each well, but in a more 
general case they are on a double tangent line (see [35]). 

The spinodal decomposition is represented in the first two graphs of Figure 2 or Figure 3. 

In longer time, the separated regions evolve to reduce their interfacial energies. These diffuse 
interfaces are shortened in an effect resembling the surface tension on a sharp interface, as the 
material fronts move to reduce their own curvature (see [11] and [31]). Finally, the solution reaches 
an equilibrium the location and form of which depend on the total initial concentration (see ]24]). 
Nevertheless this equilibrium is always a solution with an interface with minimal measure. On 
Figures 2 and 3, this phenomenon is observed on the last four graphs. 

Figure 2 corresponds to an evolution under the classic quartic double-well potential (0.9) with 
non scaled coefficients, whereas Figure 3 corresponds to an evolution under the logarithmic potential 
(0.7). They are both simulated on a 12 x 12 mesh under Qi polynomial elements. The e parameter 
is such that = 0.07. We see that the evolutions are quite similar and lead to the same stationary 
state. On these two evolutions, we can compare the difference of the energies or the norm of the 
difference (see next paragraph). Note that the polynomial approximation of the logarithm does 
not change the qualitative behavior. The same patterns appear and the long time behavior is very 
similar. The only notable difference is that with the logarithmic nonlinearity, the dynamic is slower. 
This is particularly clear on the graphs (b). The spinodal decomposition is almost completed only 
for the polynomial. Similarly, on graphs (f), we see that at time t = 1, the logarithmic evolution 
has not reached equilibrium yet. We have observed this in all our tests. 

The second evolution is often illustrated by the classical benchmark cross. It can be considered 
as a qualitative validation of the numerical methods. This long time behavior is illustrated in Fig- 
ure 4. Starting from a cross-shaped initial condition, the interface first diffuses from the arbitrary 
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Figure 1: Free energy density for two different temperatures. 



width of the initial condition to the equihbrium interface width. Next, the solution tries to reduce 
its interfacial energy and tends to a circular form. In the total free energy (0.3), the term with the 
free energy function / is responsible to the spinodal decomposition, whereas the gradient term is 
responsible for the interfacial reduction. This phenomenon has been simulated on a 256 x 256 mesh 
under Q3 polynomial elements. Figures 4 (a), (b) are obtained with the quartic nonlinearity. We 
see on Figure 4 (c) and (d) that again the qualitative behavior is very similar with the logarithm. 

It is difficult to measure precisely the qualitative difference between the two evolutions. The 
only physical quantity which can be measured in two dimensions is the energy. A detailed study 
of this aspect is performed below. Moreover in the one-dimensional case, we are able to measure 
the interface. We will see that the quartic nonlinearity tends to thicken the interface. 

The replacement of the logarithmic free energy by the quartic one has been done by many 
authors in order to avoid numerical and theoretical difhculties raised by the singular values ±1. 
More generally, we can discuss the approximation of the logarithm by polynomial functions. We 
consider the 2n-th order polynomial Taylor expansion f2n- 



hn ■.^u^\TA ^-^ ) + T 



^2p(2p-l)_ 



K2u- (2.2) 



It is defined up to an additive constant K2n- The constant if 2,1 is apparently arbitrary. However, 
it is preferable to choose it in order that the energy of a solution u 

^2„(w) / (/2„(u) + n\Vu\'') AV (2.3) 

is well defined on unbounded domains. Since it is expected that the solution converges to one of 
the binodal values, it is natural to choose so that /2„ vanishes at those points. We always 
consider this choice. 

We have seen above that the quartic approximation does not seem to change drastically the 
qualitative behaviour, except that the evolution is faster. We now perform a quantitative study to 
measure more precisely the effect of the polynomial approximation. 
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(e) t=0.6 (f) t=l 

Figure 2: Spinodal decomposition under the classic quartic double- well potential. 
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(e) t=0.6 (f) t=l 

Figure 3: Spinodal decomposition under a logarithmic potential. 
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(a) Quartic potential - t=0 



(c) Logarithmic potential - t=0 



(b) t=l 



o 



(d) t=l 



Figure 4: Evolution of a cross-shaped initial condition to a bubble. 
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The spinodal and binodal points are drawn in Figure 5 for various n. When n increases, the 
spinodal and binodal points converge to the corresponding values for the logarithmic potential. 
However, the convergence is rather slow (see Figures 5 and 6). 




— * — Polynomial binodal points 
— e — Polynomial spinodal points 

Logarithmic binodal point (i 

Logarithmic spinodal point 0.5 _ 
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Figure 5: Polynomial and logarithmic spinodal and binodal points. 

In the one-dimensional case, it is possible to study the thickness of the interface. Let us consider 
the domain 17 = R and the quartic potential 



(2.4) 



which is the derivative of 



/4 := u Tc 



T 
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(2.5) 



Then a stationary solution of the Cahn-Hilliard equation (0.6) can be explicitly computed (under 
a constant mobility M{u) = 1), see [Ki]: 



u : X ^ M+ tanh (x/i) , 



where 



^3(|-l) and 1^=^^^- 



(2.6) 



(2.7) 



It is important to remark that the solution is constrained in [— We can define a char- 
acteristic length i (see Figure 7), corresponding to the width of the region containing the main 
variations of a solution u : 



I lim;^^+oo u{x)\ + \ limx^^oo u{x) \ 
Slope in interface point 
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Polynomial degree 



Figure 6: Rate of convergence of the polynomial points. 



where the interface point is the point xq where u(xo) = 0. Thus we can compute explicitly this 
length and obtain: 

2u+ 2eV2 
2e\/2 

Cahn and Hilliard have defined a parameter A := in order to characterize the interface 

length. With this parameter A we obtain the following expression for £: 

' (2.9) 



1 - ^ 



Cahn and Hilliard have shown that in the case of the logarithmic free density the interface length 
is of the same order. This suggests that the quartic double well approximation preserves important 
features of the solution. 

In Figure 8, we present the numerical solution for fl = [0, 1] (blue stars), and the "tanh-profile" 
whose coefhcients u+ and /i have been fitted to the data. The fitting on u+ corresponds to the 
value of the solution on the boundaries of the domain f2. And the fitting on /.i corresponds to a 
least square method between the numerical solution and a "tanh-profile" solution interpolated on 
the same meshes. The "tanh-profile" (defined over R) may be considered as a good approximation 
of the solution on D, = [0, 1] since the interface is very thin. The numerical solution is computed 
with 35 Qa-elements. 

However, we have measured numerically the interface width in the quartic and logarithmic 
cases. This width is plotted for various e on Figure 9. We see that as expected by the formula 
(2.8), it varies linearly with e. But, for s not too small, the interface width is thinner for the 
logarithmic equation. The quartic approximation introduces a non negligible extra diffusivity. 

We can also compare the total free energies. Denote by u the solution of a simulation with 
the logarithmic function / and by (u2n)n>2 the family of solutions of the simulations with the 
polynomial functions (/2n)n>2- For the energy, we take as reference the logarithmic total free 
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Figure 7: Interface length for the solution U4. 



Fitted curve of the numerical solution 
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Figure 8: Fitted curve on the "tanh-profile". 
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Figure 9: Length of the interface for the quartic and logarithmic potentials, 
energy, and we study 

\F{u2n)-J'{u)\. (2.10) 

On Figure 10, the evolution of the logarithm of this quantity is plotted during a classical 
spinodal decomposition in dimension one. 

We can see important peaks at the begining and smoother peaks between iterations 500 and 
700. These peaks appear when the solution has a rapid evolution and when its topological form 
changes. For instance, these peaks correspond to the changes beetween the fourth and the fifth 
images of Figure 2, and between the fifth and the sixth images. After the iteration 750, all the 
solutions are in an asymptotic stable state, and the energies do not change anymore. 

For a quartic potential [n — 2 i.e. in Figure 10), the energy error is significant and the 
polynomial approximation is not good in that respect. 

We could as well have shown the evolution of 

\^2n{u2n) - ^ {u)\ 

In fact, it is very similar and does not bring new information. 

On a mathematical point of view, it is interesting to study the error in LF' norm: 

\u2n - u\^dx 

We see on Figure 11 that for n = 2, the error is important. It decreases with n but is still significant 
for n = 3. For n > 6, it is negligible. 

Figures 12 and 13 present the same quantities for a two-dimensional spinodal decomposition. 
We observe the same quantitative difference. Note that we clearly see that the energy evolution 
slows down as the degree n grows. 

We conclude that the classical quartic approximation of the free energy may be considered as 
a good approximation for qualitative behaviour but it produces a significant error and accelerates 
the dynamics. If precision is required, one should consider an approximation with a higher order 
polynomial. 
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Logarithmic energies under polynomial non-linearity. 
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Figure 10: Polynomial energies versus logarithmic energy. 
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ure 11: errors between the polynomial solutions and the logarithmic solution. 



14 



Two dimensional case 
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Figure 12: Polynomial energies versus logarithmic energy. 



Two dimensional case 
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Figure 13: errors between the polynomial solutions and the logarithmic solution. 
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3 Validation of the numerical method. Choice of the degree 
of the elements 



On Figure 14. we have drawn a numerical solution for different times. It is a Qi solution on a 
mesh with 100 elements under the quartic double- well potential. On figure 14(f), the solution has 
reached its stable state and has binodal values ±1 on the boundary. 




(d) t=10 (e) t=20 (f) t=50 



Figure 14: Qi solution on a mesh with 100 elements. 

In our first set of tests, we start with the same initial state near the "tanh profile" solution. 
The evolutions are driven by the quartic potential function. We wait for the stabilization of all the 
solutions and study the error on the energies and on the slopes of the interface. 

Remark that we can explicitly compute the energy of the explicit solution. And since the energy 
of a numerical simulation is decreasing in time, this energy should converge to the energy of the 
explicit solution. Figure 15 shows the evolutions of the errors between the numerical energies and 
the explicit energy according to the degree of the polynomial space P. Before the 800th iteration in 
time, the solutions are not stable. They try to minimize their energies. After the 800th iteration, 
all the solutions are in a stable state. We can see that the evolutions are qualitatively similar at 
the beginning, but the elements Qi, Q2 and don't achieve the tolerance zone, whereas the other 
elements do. However, Q2 and give a very good result. 

The slope of the interface is an essential physical quantity. So we have compared the errors 
on the slopes between the numerical solutions and the theoretical solution. Note that these slopes 
correspond to the values of the derivatives of the numerical solutions at the interface and our finite 
elements have not a C ^ regularity. 

Under the quartic double- well potential (2.5), we have an explicit slope /i for the stationary 
solution. On Figure 16, we present the numerical solution for fl = [0, 1] (blue stars), and the "tanh- 
profile" whose coefficients and /i have been fitted to the data. The fitting on corresponds to 
the value of the solution on the boundaries of the domain CI. And the fitting on /i corresponds to a 
least square method between the numerical solution and a "tanh-profile" solution interpolated on 
the same meshes. The "tanh-profile" (defined over M) may be considered as a good approximation 
of the solution on 51 = [0, 1] since the interface is very thin. 

If we want to compare the solutions between a Qi simulation and a Qio simulation, we need 
to compare the two simulations under a same complexity which, up to the inversions of the linear 
systems, corresponds to a similar computational cost. In the one dimensional case, the complexity 
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Decimal logarithm of the error on the energy. 
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Figure 15: Energies during an evolution. 

corresponds to the value Degree x Number of elements. For a Qio simulation, we only need a 
mesh with 10 times less elements than for a Qi simulation. 

Figure 16 represents the numerical solution over mesh grids with three different complexities 18, 
36 and 72, and under polynomial functions of degree 1, 2 and 3. For instance, for the elements (52, 
it corresponds to the mesh grids with 9, 18 and 36 elements. If we increase the number of elements 
or the degree of the polynomial space P, then we obtain a better approximation of the slope of the 
"tanh-profil" solution. But for the same complexity, the curves are qualitatively similar. Figures 
17(a) and 17(b) show the evolution of this approximation error according to the complexity for Qi, 
Q2 and Q3 simulations. On Figure 17(b), we have used a logarithmic scale in order to compare 
the rate of the convergence. 
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(a) Mesh 18 - Ql 



(b) Mesh 36 - Ql 



(c) Mesh 72 - Ql 




(g) Mesh 6 - Q3 (h) Mesh 12 - Q3 (i) Mesh 24 - Q3 

Figure 16: Fitted curves on the "tanh-profile". 
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(a) Versus the number of elements of the mesh (b) Versus the logarithm of the number of ele- 

ments of the mesh 

Figure 17: Comparison between the errors on the slope under the same complexity. 
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We obviously conclude that, for elements Qi, Q2 or Q^, a fine mesh allows a better approxima- 
tion. But the Q2 and Qs elements seem to reach faster a saturation. They only need 500 elements 
in order to reach a 10~^ precision, whereas the Qi elements need 5000 elements ! Figure 17(b) 
highlights this better speed on the approximation error of the slope. But Q2 and elements 
seems to have a similar speed before reaching the saturation zone. 

If we fix the complexity, we can test which degree of the polynomial space P can provide the 
best speed. Figure 18 shows this approximation error according to the degree of the polynomial 
space P under the same complexity - quantified by the number of degrees of freedom (DoF) of the 
finite elements space. 




(c) DoF 630 



Figure 18: Error on the slope versus the degree of the polynomial space P under the same com- 
plexity. 

Under the same complexity, we see on Figure 18 that high degree elements still provide better 
approximations than Qi elements. Although very high degree elements always provide better 
approximations than low degree elements, the slopes on Figure 18(c) of the curves for low degrees 
suggest that Q3 elements are a good choice. Higher elements increase the computation time for 
matrix inversion and the gain is not valuable. 

Figure 19(a) shows the error on the energies according to the complexity under Qi, Q2, Q3 and 

elements. As for the slopes, we see that the error is decreasing as the number of elements of 
the mesh is increasing. Whereas the error reaches a 10~* precision for the slopes before saturation, 
the error on the energy reaches the tolerance zone for (54 elements on a mesh with 500 elements. 
Figure 19(b) shows the logarithm of the error according to the logarithm of the complexity. We 
see that the evolution is linear for the finest meshes with a good speed. We conclude in particular 
that we can compute an order of the speed of the convergence. For Qi elements, we find an order 
2, for Q2 elements, we find an order 4, for elements, we find an order 4 and for (54 elements, 
we find an order 6. Note that the error on the energies should be of the order as the error. 

Now, we fix the complexity and compare the approximation error on the energy according to 
the degree of the polynomial space P. For the complexities 90, 180, 300 and 630, we have drawn 
the decimal logarithm of the errors on Figure 20. 

Again, under a same complexity, if we increase the degree of the polynomial space P, the high 
degrees can provide better approximation, except on the coarse grids. We can conclude that for a 
fixed mesh (fine enough), high degrees provide a better approximation. But for each complexity, 
it seems that we have a saturation because the Qe, Qi, Qsi Q9 and Qio elements have almost the 
same errors. We conclude that we have to use elements with high degrees, but it is not necessary 
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TalerancB of Newton Iterates 
Tolerance System Inversion 

200 400 600 800 1000 1200 1400 1600 1800 2000 



Tolerance of Newton I 



2.4 2.6 2.8 



(b) Errors according to the logarithm of the 
complexity 



(a) Errors according to the complexity 

Figure 19: Errors according to the logarithm of the complexity 
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to choose the highest. We have to take into account the computational cost, and the precision of 
our inverse solver. Indeed, even if the complexity is the same, the finite elements matrices have not 
the same profil. For instance, the bandwidth of the "mass" matrix for Qio elements is much larger 
than for Qi elements. Figures 19 and 20 indicate that Q2 and Q3 elements are a good compromise 
to ensure good results without increasing the computational cost too much. 

In the two dimensional case, the results are drawn on Figure 21. The behaviour is similar. 



-2.5 r 




xlO* 3.2 3.4 3.6 3.8 4 4.2 4.4 4.6 4.8 5 



(a) Errors according to the complexity (b) Errors according to the logarithm of the 

complexity 



Figure 21: Errors on the energy according to the logarithm of the complexity 

The energy and the interface are essential physical quantities. From a mathematical point of 
view, it is also important to study the error. 

Figures 22(a) and 22(b) show the error according to the complexity for Qi, Q2, Qa, Qi 
and Q5 elements. On Figure 22(b), we have used a logarithmic scale in order to compare the 
convergence rate. We have computed the order of the speed of the convergence. If we extrapolate 





200 400 600 800 1000 1200 1400 1600 1800 2000 

(a) Versus the complexity. 



(b) Versus the decimal logarithm of the com- 
plexity. 



Figure 22: Comparison between the errors under the same complexity, 
the lines, we can find the necessary complexity in order to reach the saturation. 



Degrees Order Complexity for saturation Grid for saturation 

1 1.9960 423829 423829 

2 3.9682 3558 1779 

3 4.0216 4110 1370 

4 4.9119 2364 591 

5 5.9040 1475 295 



Again, Q2 and elements give very good results for a reasonable computational cost. We have 
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decided to prefer elements because it seems that they provide better results on the interface 
length as shown on Figure 18. 



4 Stationary states 

The Cahn-Hilliard equation has a lot of asymptotic equilibria (see [18], [29] and [30]). In the one 
dimensional case, a state can be described by the number of interfaces and their positions. On 
Figure 23, we show four states which are numerically stable. It is possible to observe more than 
one interface only for small e. Only when the interface is very thin - i.e. for small e, the interfaces 
do not interact. Note that the energy increases with the number of the interfaces. 

In fact, this is a bifurcation phenomenon. When e crosses critical values, bifurcations happen 
and more stationary solutions appear. 



One interface 



0.2 0.4 0.6 0.8 



Two interfaces 



0.2 0.4 0.6 0.8 1 



Three interfaces 



0.2 0.4 0.6 0.8 f 



Four interfaces 



0.2 0.4 0.6 0.8 1 



Figure 23: Four numerically stable states. 



In [24[, the authors consider the stationary states of (0.6) on the square. They numerically 
study the solutions of the following semi-linear elliptic equation. 



c = u — + e^Au, on fi, 



Vu • = 0, 



on 951, 



(4.1) 



together with the mass constraint: 



^/^u(.)dx = ^, 



(4.2) 



where = [0, 1]^ is the square, c G M and to e M are parameters. They study stationary solutions 
under the three-dimensional parameter space (c, to, 1/e^). For this system and for all e, a trivial 
solution is given by the constant solution u = m with c = m — vn? . The linearization around u = m 
of (4.1) under the mass constraint reads 



= (l - 3?7i2) u + e'^Au, on fl, 
Vu ■ v = 0, on (917. 



(4.3) 



Let Vr be an eigenfunction of the Laplace-Neumann operator in 17 defined in (4.6) with eigenvalue 
r E M^, then Vr is also an eigenfunction of (4.3) when 



1 



1 -3to2 



for |7Tl| < 



^/3■ 



(4.4) 
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But (T_|_ = 1/ for the quartic double-well potential, so this equality shows that bifurcations may 
occur only for m in the spinodal region. For the square domain = [0, 1]^, the eigenfunctions are: 

Vr{x,y) =Vk,i{x,y) := cos^irkx) cos{T:ly) for {x,y) G [0,1]^, (4.5) 

with (fc, /) £ such that r = (fc^ -I- 1^)tt'^ . For the mode vi^i (i.e. r ~ 27r^), we obtain non- 
trivial solutions bifurcating at u = ±m* with m* = ^^(1 — e^r) /3. We fix m = 0, such that the 
bifurcations occur as 1 = e^r. 

The previous asymptotic equilibria - described in [24] - are asymptotic solutions of the dynam- 
ical evolution. For instance, we have obtained the vi^i mode as a stationary solution of a dynamical 
evolution (See Figure 24(c)). A random start may lead to different modes, and actually we only 
see the most stable of them in long time. Figures 24(a) and 24(b) show the stable states that we 
see most of the time. 

All the symmetrical states are also stable. In [23], the authors have studied the global attractor 
on a square and they have proved that, after the first bifurcation, there exist 4 minimal attractors 
(see Theorem 4.2 in ]2.3]) obtained by symmetrization of Figure 24(a). The other stable states 
shown here appear after subsequent bifurcations. Starting the simulation with well chosen initial 
data, we have been able to recover dynamically all the stable states described in [24]. If we 
choose the mode V4^i +vi^4 (which is the last mode studied by Maier-Paape and Miller), we see on 
Figure 24(d) the asymptotic equilibria that we have obtained. 
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In general, the stable states are deeply dependent on the eigenvalues of the Laplace operator 
on il. Let {pk)^.^f>j and (wp^)j,gf^ be the eigenvalues and eigenvectors of the following problem: 



( -^^'Pk = PkVp^, on n C 
Vwpj^ ■ ly — 0, on dfl, 



(4.6) 



In [23], the authors have studied the bifurcations and the global attractors of the Cahn-Hilliard 
problem. In their nomenclature, they consider the following Cahn-Hilliard equation: 



f dtv = Aw, on SI CM", 

w = -Xv + 72W^ + 73W^ - Av, on C M", 
Wv • = = • i^, on dft, 



(4.7) 



where A, 72 and 73 are parameters. If u is a solution of the system (0.6) on f2, then w is a solution 
on Sl/Ve of (4.7) if we define for alH G M and x e fi: 



V : {t,x) i-> u{t,x\/e). 
and the correpondence is given by the following equalities. 

A := -, 72 := and 73 := -. 



(4.8) 



(4.9) 



They prove that the first bifurcation occurs as their parameter A is greater than a particular value. 
For our problem, this bifurcation occurs as — > pi. 

Below, we study this first bifurcation and illustrate theoretical results of [23]. 



4.1 Asymptotic stable states on a rectangle 

In the case of a rectangular domain H = [0, 2] x [0, 1], the hypothesis of the Theorem 4.1 in [23[ 
1 7r2 

holds. Accordingly, if — > pi := — then there exist exactly two attractors which can be 
E'^ 4 

expressed as 



2£ / 1 7r2 



± We (a;, y) = ±^ Y ^ - — cos 



2 / 



) +y/e o 



1 



1/2N 



xe[0,2],ye[0,l]. (4.10) 



We define the approximated attractors ±11^ by 



1 TT 



±v,{x) :=±C(£)Y^-— cos(^— j, xe [0,2], ye [0,1], 

where C{e) is a constant depending on e. It is chosen in order to minimize the L^ norm of u'^ — v^. 

2 

For multiple values of the parameter e around the value — , we have obtained the corresponding 

TT 

numerical stationnary states u'^ . 

We have checked numerically the validity of formula (4.10). We study the following quantity 



\\Veh 



This relative L^ norm should converge to zero. On Figure 25, we have plotted the decimal logarithm 

2 . . . 1 . . . 

of this relative L norm according to the decimal logarithm of — . Figure 25 is in conformity 
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Figure 25: Convergence of the "bifurcationned" solutions on a rectangular domain. 




Figure 26: Order of convergence of the minimal constant C(e). 



„ 1 

with the expecting theoretical results. We can see that the relative L error converges to as — 



converges to — . We even can improve formula (4.10) and find the exponent arectanguiar such that 



TT 



arectangida,- 



where C is an unknown constant. We find that the exponent arectanguiar — 1/2 + 0.98908, almost 
3/2. Moreover, since we know explicitly the attractor, we can verify that our minimal constant C{e) 
2e 

is near — r=. On Figure 26, we have drawn the logarithm of our minimal constant C(e) according 
v3 

to the logarithm of e. 
We find 

C(e) - 1.0682 * £0-83603 (4 ^^1) 



this is in conformity with the fact that C(e) converges to lim 



2e 



The segment may be seen as a degenerate rectangle. On the segment [0, 1], the eigenvalues of 
(4.6) are pk := fc^Tr^ for all fc e N. According to Theorem 4.2 in [23], there is a bifurcation at 
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— > TT . Moreover. Remark 4.2 in states that there exist two minimal attractors ±Mg which 
can be expressed as 



± Us{x) = ±C(e)i / ( -J ^ ""^ I cos (ttx) + Ve o 



1/2N 



e[0,l], (4.12) 



where C(£) is a constant which can depend on e. Again, we define the approximated attractors 
±We by 

±Ve{x) := ±C{e)J (\~nA cos (ttx) , x € [0, 1]- 



For multiple values of the parameter e around the value — , we have obtained the corresponding 

TT 

numerical stationary states u'^. We choose the constant C{e) in order to minimize the norm of 
— and study the convergence of u'^ to using the quantity 

||<-«,!|2_ 2||<-t;,||2 



Figure 27: Convergence of the "bifurcationned" solutions on a segment. 



Figure 27 is in conformity with the expected theoretical results. We can see that the relative 
1 



L error converges to as — converges to tt and find the exponent a 



segment 



such that 



where C is an unknown constant. We have found asegment = 1/2 + 1.0254, again almost 3/2. 
We have said that the constants C(e) have been numerically chosen in order to minimize the 

2 



L norm of u'^ — v^. If we extend the results of [2.')], we expect that the constant C(e) 



Thus, on Figure 28, we have drawn the logarithm of our minimal constant C{e) according to the 
logarithm of e. If we study the slope, we find 



C{e) - 1.0844 * £"-9459* 
which is again in conformity with the theoretical formula. 



(4.13) 
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-0.66 -0.64 -0.62 -0.6 -0.58 -0.56 -0.54 -0.52 

Figure 28: Order of convergence of the minimal contant C{e). 
4.2 Asymptotic stable states on smooth domains 

For a smooth domain il, the hypothesis of Theorem 3.1 in [23] holds. We have considered an 
ellipse. On the ellipse, the first eigenvalue pi ~ 0.8776 is simple. On Figure 29, we have drawn 

the corresponding first eigenvector. If — > Pi, the problem (4.7) has two steady states ztu^ which 



Figure 29: First eigenvector on the ellipse. 



can be expressed as 



zbuj = ±C{s)J i ^ - Pi]vp^ + y/e o 



1 



1/2N 



where C(e) is a constant which can depend on s. We define the approximated attractors ztug by 



±v^ := ±C{e)J 3 - Pi fpi, 



27 



where Vp^ is a fixed eigenvector. On the Figure 30. we have drawn a steady states. For muhiple 



5.0000E-01 

1 .0000 E+00 




Figure 30: Steady state on the ellipse. 



values of the parameter e around the value , we have obtained the corresponding numerical 

statiomiary states u'^. As in section 4.1. we choose the constant C(e) in order to minimize the 
norm of u'^ — and study the convergence of u'^ to v^. We consider the quantity 




According to theorem 3.1, this relative norm should converge to zero. On Figure 31, we have 

2 ■ ■ ■ 1 

drawn the decimal logarithm of this relative L norm according to the decimal logarithm of — — pi. 

-1.2r 




Figure 31: Convergence of the "bifurcationned" solutions on the ellipse. 
Figure 31 corroborates the expected theoretical results. We can see that the relative error 
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Figure 32: Eigenvector of the first eigenvalue on the trapezoid. 



converges to as — converges to pi . Then we compute the exponent aeiupse such that 



\\u',-v,h^C 



where C is an unknown constant. We find that exponent aeiupse = 1/2 + 0.9580. 



4.3 Asymptotic stable states on a trapezoid 

We try to see if the results of [2.!] extend to non smooth domains. We have tested a trapezoid 
where the first eigenvalue pi ~ 2.2417 is simple. On Figure 32, we have drawn the corresponding 
first eigenvector. 

The two steady states ±.u^ should be expressed as 



Pi Vpi + V £ o 



72 



1/2N 



x e [0,1], 



where C{e) is a constant which can depend on s. We define the approximated attractors ztw^ by 



1 



- Pi ]vpi, X e [0,1] 



where Vp-^ is a fixed eigenvector. On Figure 33, we have drawn a steady state. 

For multiple values of the parameter e around the value , we have obtained the correspond- 

VPi 

ing numerical stationnary states u'^. As in section 4.1, if we choose the constant C(e) in order to 
minimize the norm of u'^ ~ v^, we can study the convergence of u'^ to v^. We have found that 



IK - Veh 



does not converge to 0. It seems that the bifurcation is different in this case. On Figure 34, we 
have drawn the decimal logarithm of \\u'^ — Ve\\2 according to the decimal logarithm of — — pi We 
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find that 



atrapczoid 



Pi 



with atrapezoid = 0.49937. Tlius this difference is of the same order as each term. 

As in the case of the square, we can find numerically stable states corresponding to the next 
modes in the nomenclature of Maier-Paape and Miller in [24] . We have found 4 numerically stable 
states (see Figure 35), with energies that have been drawn on Figure 36 against the length of 
the interface. We can clear see the linear dependance between the two (the red line is the linear 
regression according to the least square method). 




Figure 35: Numerically stable states on the trapezoid. 
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Energies according to the interface length 

-0.53 

-0.54 

-0.55 

-0.56 

-0.57 
at 

m -0.58 
-0.59 
-0.6 
-0.61 



"■"1.5 1.6 1.7 1.8 1.9 2 2.1 2.2 2.3 

Interface length 

Figure 36: Energies of the numerically stable states on the trapezoid according to the length of 
the interface. 
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